Repository navigation
ForceFreeStates.Galerkin - FEATURE! - Solve the Galerkin system by banded Cholesky, with LU fallback - #513
Conversation
…r gal_solver = "cholesky"
gal_solver = "cholesky" called LinearAlgebra.LAPACK.pbtrf!/pbtrs!, which do not exist, so the
option died with an UndefVarError. Call LAPACK zpbtrf('L')/zpbtrs('L') directly through
libblastrampoline, unscaled, on the lower band the assembly already builds for Cholesky. Unlike
the reference implementation, a non-positive-definite matrix (zpbtrf info > 0) raises an error
naming gal_solver = "LU" instead of continuing with a partial factor.
Add examples/LAR_gal_cholesky_example (LAR analytic equilibrium, positive-definite Galerkin
matrix), regression case lar_gal_cholesky, test/runtests_galerkin_cholesky.jl, and a solver note
in docs/src/galerkin.md.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
matt-pharr
left a comment
matt-pharr
left a comment
There was a problem hiding this comment.
Couple cleanup items and one main question, being whether we should do away with the flag altogether and just make every galerkin solve try to use cholesky first and then use LU as a fallback. See the relevant comment.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…to enable Cholesky solver for ncoil > 0
…rst and fall back to LU gal_solver defaults to "auto", following @matt-pharr's review recommendation on #513: the full-band Galerkin matrix is solved by banded Hermitian Cholesky when it is Hermitian positive definite (Dirichlet rows' columns moved to the right-hand side), else by the Jacobi-scaled LU. "LU" skips the attempt; "cholesky" errors instead of falling back. The matrix is always assembled full band, so the lower-band assembly branches are gone. The test now checks the solver choice and the solution against a dense solve. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…try by entry, diagonal included The Hermiticity check before the banded Cholesky skipped the diagonal and used one tolerance, 1e-12 of the largest entry. zpbtrf reads only the real part of the diagonal, so a complex diagonal (the resonant elements' conj(u)·Lu) was dropped silently; and since the diagonal spans many decades, the small rows passed any asymmetry. Each entry is now compared with 1e-10·√|aᵢᵢaⱼⱼ|, diagonal included (the quadrature round-off is ~1e-11 of that scale). gal_cholesky_solve! requires kl == ku and solves in place of sol once zpbtrf succeeds. Galerkin() now defaults to solver = "auto" like gal_solver, and its docstring lists "auto". Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…holesky falls back to LU only when not positive definite Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
StuartBenjamin
changed the title
github-actions
Bot
added
the
changed-results
…sonant surfaces Resonant powers are the eigenvalues of each side's own M₀ block. √(−det) matches those eigenvalues only when the trace vanishes, and that path is unchanged. alpha_override and extended_precision_bvp are removed. Co-authored-by: Cursor <cursoragent@cursor.com>
… bitwise-regression for simplicity and correctness
…sonant surfaces (#517)
d-burg
left a comment
d-burg
left a comment
There was a problem hiding this comment.
#517 was merged into this branch, so I'm reviewing its exponent change here, where it will reach develop.
The exponent change is correct, and it is the paper's definition rather than a departure from it, so that's good. Glasser (2016) defines the powers by det[M₀ − αI] = 0 (Eq. 46). Of α = ±√(−D_I) (Eqs. 48–49) the paper says it could not be proved analytically, only that the roots converge to it as Fourier components are added, and it proposes that agreement as a Fourier-convergence check. So sqrt(−det) was Eq. 46 under an unstated tr = 0 assumption.
From claude: regarding @matt-pharr's question (more right, or only more self-consistent): as far as I can derive it, the resonant block has the form [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real, so tr = −2i·Im κ. That is a genuine term of the truncated operator, not an assembly slip. Matching to a basis with any other powers makes Δ′ depend on the matching distance, so this is more right for the truncated system. It does hide one thing: the old imaginary Δ′ was the only visible symptom of a non-converged resonant block, so I'd like the trace itself reported. Log Im tr / 2 and α per surface (ideally also next to alpha in the HDF5 output), and warn when |tr| is not small against α. The debug print in Asymptotics.jl (lines 125–126, outside this diff) still prints di next to an α that is no longer √(−di); something like:
block_det = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2]
block_tr = m0mat[1, 1] + m0mat[2, 2]
msg *= @sprintf(" det= %+.12e, tr= %+.12e %+.12ei\n", real(block_det), real(block_tr), imag(block_tr))
msg *= @sprintf(" alpha= %+.12e %+.12ei\n", real(alpha[1]), imag(alpha[1]))Requests before this goes to develop:
- Release note. This PR now carries #517's results-moving change and the removal of
extended_precision_bvp, but the release note covers onlygal_solver. Suggested Migration: "To reproduce the previous Galerkin results, setgal_solver = "LU". Removeextended_precision_bvpfrom[ForceFreeStates](the TOML key is ignored with a warning; passing it as a keyword tosolveis now an error); the Δ′ BVP always solves in Float64. Anything pinned to imaginary Δ′ or to Galerkinalphamust be re-baselined." The Numerical-impact line needs the exponent change too. - CI is red at
test/runtests_eulerlagrange.jl:706: the opposed-powers testset still checks bitwise a branch that was removed (suggestion inline). - Harness at the current head. The stamp here (
cd290c659) predates the #517 merge, and #517's own report was against a merge-base older than the axis-start change (#514). The near-axis surface is the one most exposed to that. - Wording: "Hermitian with a real diagonal" rather than "Δ′ must be real", as Matt suggested on #517.
- Solovev. Solovev is up-down symmetric, so its blocks should be traceless and the exponent change should not move it, yet #517 reported Im Δ′ going from −1.8e9 to −27 there. What is Im tr at its q = 2 surface? If it is round-off level, that change needs another explanation; if not, the matrices have an imaginary part on a symmetric equilibrium.
- Near-axis surface at n = 4. A residual Im Δ′ of −25 (about 1 % of the diagonal) means the series and the integrated solution disagree by about 1 % there; I'd treat that as the error bar. A cheap diagnostic: run with
JULIA_DEBUG=GeneralizedPerturbedEquilibrium, compare Im tr, α and dpsi^(−2α) at q = 1.25 with the other surfaces, then a three-pointsingfac_minscan (3e-4, 1e-4, 3e-5) to see whether Δ′₁₁ plateaus. - Small:
using DoubleFloatsinsrc/ForceFreeStates/ForceFreeStates.jl(line 22) is now dead; GGJ has its ownusing.
| @testset "column scaling stays /dpsi^α when the powers are opposed" begin | ||
| α = ComplexF64(1.5) | ||
| dpsi = 1.3e-4 | ||
| big = ones(ComplexF64, 2, 2) | ||
| small = ones(ComplexF64, 2, 2) | ||
| FFS._apply_column_powers!(big, small, -α, α, dpsi) | ||
| pfac = dpsi^α | ||
| big_ref = ones(ComplexF64, 2, 2) | ||
| small_ref = ones(ComplexF64, 2, 2) | ||
| big_ref ./= pfac | ||
| small_ref .*= pfac | ||
| @test big == big_ref | ||
| @test small == small_ref | ||
|
|
There was a problem hiding this comment.
The opposed-powers branch this testset pinned is gone, so its bitwise check now fails CI (line 706). The second half of the testset already checks the general path exactly, so I'd drop the first half rather than relax it to ≈.
| @testset "column scaling stays /dpsi^α when the powers are opposed" begin | |
| α = ComplexF64(1.5) | |
| dpsi = 1.3e-4 | |
| big = ones(ComplexF64, 2, 2) | |
| small = ones(ComplexF64, 2, 2) | |
| FFS._apply_column_powers!(big, small, -α, α, dpsi) | |
| pfac = dpsi^α | |
| big_ref = ones(ComplexF64, 2, 2) | |
| small_ref = ones(ComplexF64, 2, 2) | |
| big_ref ./= pfac | |
| small_ref .*= pfac | |
| @test big == big_ref | |
| @test small == small_ref | |
| @testset "column scaling multiplies each column by dpsi^power" begin | |
| dpsi = 1.3e-4 |
There was a problem hiding this comment.
Done in a968945: the opposed-powers half is dropped, and the general half stays.
There was a problem hiding this comment.
Confirmed in the final head: Stuart removed the obsolete opposed-powers half, while retaining the general column-scaling assertions. We kept that change through the develop merges. See the current test. The final post-merge regression run also completed successfully: 21 cases at three refs, all 63 runs.
| function resonant_block_exponents(m0::AbstractMatrix, sig::Real) | ||
| tr = m0[1, 1] + m0[2, 2] | ||
| det = m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2] | ||
| α = sqrt(tr * tr / 4 - det) |
There was a problem hiding this comment.
With the traceless branch removed this is the only path, so a real block with complex roots throws DomainError and one with real roots returns Float64. Converting keeps it total and type-stable without changing any complex result.
| α = sqrt(tr * tr / 4 - det) | |
| α = sqrt(ComplexF64(tr * tr / 4 - det)) |
There was a problem hiding this comment.
Done: α = sqrt(ComplexF64(tr * tr / 4 - det)).
There was a problem hiding this comment.
Retained Stuart’s fix: resonant_block_exponents computes sqrt(ComplexF64(tr * tr / 4 - det)), including for real inputs with complex roots. See the helper.
| function resonant_block_exponents(m0::AbstractMatrix, sig::Real) | ||
| tr = m0[1, 1] + m0[2, 2] | ||
| det = m0[1, 1] * m0[2, 2] - m0[2, 1] * m0[1, 2] | ||
| α = sqrt(tr * tr / 4 - det) |
There was a problem hiding this comment.
When tr²/4 − det is negative (a Mercier-unstable surface), it is real by construction, so the sign of the root comes from the sign of a zero or round-off imaginary part. Left and right can then label opposite solutions "big"; a standalone copy of this helper does that for a real block and its mirror. The old shared α made the pairing consistent by construction, and it is now gone entirely. Could α be canonicalised the same way on both sides (Re α > 0, or Im α > 0 when Re α is zero to round-off), or the pair computed once per surface and passed to both sides?
There was a problem hiding this comment.
Retained the canonicalization: when the real part of alpha is roundoff-sized, we choose nonnegative Im(alpha), and use the direction-aware tie-breaking for the big/small roots. The Mercier-unstable test perturbs the radicand’s imaginary sign and checks that alpha and the big-root choice remain consistent. See the tests. Each side still uses its own block’s eigenvalue powers; the old alpha_override is not restored.
| end | ||
| end | ||
|
|
||
| @testset "resonant block Frobenius exponents" begin |
There was a problem hiding this comment.
These tests pin the 2×2 algebra but not the claim. The assembled block always has a purely imaginary trace, which none of the inputs here has, and the unit-test fixtures are up-down symmetric, so CI never exercises the new path. I'd add a case with the real block structure (before the column-scaling testset):
@testset "Hermitian-structured block: shared imaginary shift, real eigenvectors, same powers on both sides" begin
# Shape the assembly produces: [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real and κ → −κ on the left side.
κ, φ, γ = 0.3 + 0.02im, 0.7, 2.5
m_right = ComplexF64[-κ+0.5 φ; γ conj(κ)-0.5]
m_left = ComplexF64[κ-0.5 φ; γ -conj(κ)+0.5]
α_right, big_right, small_right = FFS.resonant_block_exponents(m_right, 1.0)
α_left, big_left, small_left = FFS.resonant_block_exponents(m_left, -1.0)
@test α_right ≈ sqrt((real(κ) - 0.5)^2 + φ * γ)
@test α_left ≈ α_right
@test imag(big_right) ≈ -imag(κ)
@test imag(small_right) ≈ -imag(κ)
@test big_left ≈ big_right
@test small_left ≈ small_right
# ξ′ entries of both leading eigenvectors are real, so the big/small cross-Wronskian 2α/φ is real.
@test abs(imag(-(m_right[1, 1] - big_right) / φ)) < 1e-14
@test abs(imag(-(m_right[1, 1] - small_right) / φ)) < 1e-14
@test -(m_right[1, 1] - small_right) / φ + (m_right[1, 1] - big_right) / φ ≈ 2 * α_right / φ
endIdeally also one integration-level assertion on an asymmetric fixture: left and right powers agree, and the Riccati Δ′ diagonal is real to a stated tolerance.
There was a problem hiding this comment.
Added as you wrote it, before the column-scaling testset (a968945).
There was a problem hiding this comment.
Stuart added the Hermitian-structured block test you suggested, and it remains in the final head: shared imaginary shift, left/right powers, real leading eigenvector components and the cross-Wronskian check. See the test. Your additional ideal integration-level assertion on an asymmetric fixture was not added as a new strict convergence guarantee. Instead we ran the asymmetric DIII-D LU/Cholesky/Riccati comparison across resolutions and matching controls; it exposes a residual Galerkin–Riccati discrepancy rather than supporting a claim of complete agreement. That limitation and the further accuracy investigation are recorded in #603 and the PR body.
| M_solve = Tc.(M) | ||
|
|
||
| M_lu = lu(M_solve; check=false) | ||
| M_lu = lu(M; check=false) |
There was a problem hiding this comment.
I don't think the evidence supports deleting extended_precision_bvp before this reaches develop, only changing its default. It covers only singfac_min = 1e-4; #517's own 2 % miss at singfac_min = 1e-6 (Solovev q₀ = 1.85) is a case where Float64 is measurably wrong, and the code can no longer detect or repair it. A singfac_min scan is also exactly how one would verify the new exponents. The exponent fix and the precision of this LU solve are separate error sources. Could we keep the option as opt-in (default false, which keeps the whole speedup), or replace it with one step of iterative refinement (Float64 LU, residual in Double64) plus a warning when the correction to the read-off entries is large? What do you think?
There was a problem hiding this comment.
Agreed. It is back as an opt-in, default false (2c02244), so the speedup stays and a singfac_min scan can still check Float64 (point 6 is that case). A test checks that the Double64 solve matches the Float64 Δ′ diagonal on the DIII-D BVP case.
There was a problem hiding this comment.
Agreed that the exponent correction and BVP precision are separate. We retained extended_precision_bvp and restored the prior true default rather than shipping a default precision reduction here. Explicit false remains available. The no-control path also retains Double64. We added full complex-matrix and componentwise Float64/Double64 checks, plus an independent 256-bit reference for the no-control helper. No iterative-refinement replacement or default precision flip is included in this PR.
…surface; trace in the debug print; review tests - resonant_block_exponents: radicand in ComplexF64; a round-off imaginary part is dropped, so a negative tr²/4 − det gives the same α (and big root) on both sides. - Debug print: det, tr and α instead of di. - Tests: column-scaling testset keeps only its general half (the opposed-powers branch is gone); Hermitian-structured block; Mercier-unstable block with ±1e-17i in the radicand. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…impler α canonicalisation
- extended_precision_bvp (default false) solves the Δ' BVP and its PEST3 combination in
Complex{Double64}, to check Float64 at small singfac_min; it is no longer a dropped key.
- resonant_block_exponents: when Re α is round-off, Im α > 0 (same α on both sides).
- DIII-D BVP test: the Double64 solve matches the Float64 Δ' diagonal.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…t small against α Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
Thanks for the careful review. Most of it concerns #517, which is merged into this branch, so I've answered those points here too. Pushed c827d4d (three commits on top of 4b44dcd). Point by point:
My three commits on top of 4b44dcd change none of these numbers. With
You also asked for the trace to be reported. The debug print now shows det, tr and α. A warning fires when a resonant block's |tr|/2 exceeds 1e-3·|α|. On the n = 4 and Solov'ev decks above the ratio is at most 3e-5, so it stays quiet there. |
|
Just human stuart adding here: @jhalpern30 most of these comments refer to ur original PR, I tried my best to respond to d-burg's review, but feel free to check out the responses if you want |
…v; conflicts take OpenFUSIONToolkit#513's newer OpenFUSIONToolkit#517 state
|
@d-burg @StuartBenjamin Following up on point 5 and on what τ = Im tr(M₀)/2 actually is. I ran a few more tests on feyn. In short: τ is numerical, it has two identifiable sources, and both go away with resolution. The eigenvalue exponents are still the right choice. What τ is. In every run, τ equals Im[(F̄⁻¹K̄)ᵣᵣ]/(n q′) at the rational surface (r = the resonant mode) to within 1e-9. That combination is large away from the surface and should pass exactly through zero at ψₛ, because two terms cancel there: the resonant term (F̄⁻¹)ᵣᵣ·Im K̄ᵣᵣ and the sum over the other modes. τ is whatever is left over from that cancellation. It isn't geometry. Mirroring the DIII-D-like g-file in Z does not flip the sign of τ (the mirror/original ratios are +0.2 to +3.2), and the up-down symmetric Solovev has the largest τ of the cases I ran. The off-diagonal Δ′ phases, on the other hand, conjugate exactly under the mirror (e.g. Δ′₁₂ −16.57° → +16.57°). So the complex off-diagonals are genuine geometry: Δ′ is Hermitian, and it is real only for up-down symmetric equilibria. Where the leftover comes from (DIII-D-like deck):
Neither mtheta (256 to 1024) nor etol (1e-10 to 1e-14) changes either part. Solovev (point 5). Its τ (1e-6 to 1e-4) wanders and flips sign as mpsi and mtheta change, and the Mercier D_I itself drifts with the Solovev grid. So the floor comes from the analytic Solovev equilibrium rather than from the ψ/θ grids. Amplified near the pole, that gives the old −1.8e9, which fits Stuart's numbers. What it means for this PR. At any finite resolution the discretized operator really does have this small trace. The eigenvalue exponents keep the asymptotic series consistent with it, so Δ′ no longer depends on the matching distance. Forcing tr = 0 would undo that. Nothing suggests Glasser's ±√(−D_I) is a symmetric-only result; it is the converged limit. On the new warning. τ is only a loose proxy for an under-resolved surface. At mpsi = 128 the q=4 surface has det(M₀) 34 % off the Mercier D_I, yet τ/α ≈ 4e-5, below the 1e-3 threshold. Comparing det(M₀) directly with the Mercier D_I would catch that case. I'd suggest two follow-ups for an issue rather than this PR: that check, and inserting each rational surface into the radial grid as a knot, which removes the interpolation part of τ by construction. |
d-burg
left a comment
d-burg
left a comment
There was a problem hiding this comment.
Thanks @StuartBenjamin, this is a nice cleanup. Moving the Dirichlet columns to the right-hand side is the right way to keep the system Hermitian and still honor the rpec coil sources, checking zpbtrf's info is a real fix, and one full-band matrix for both paths removes a lot of branching. Before it goes in:
- Check Hermiticity before using Cholesky.
"auto"currently returns a Cholesky answer for any matrix whose lower band is positive definite, even when the full matrix isn't Hermitian, so it can silently solve a different system. Suggestion 1 adds the check (it also catches a missed Dirichlet row), and suggestion 3 is a test that fails on the current head. - Make the solver choice visible. Cholesky and LU differed by 0.40 % on the PEST3 Δ diagonal before the M₀-block fix, so a scan that crosses the positive-definite boundary would jump with nothing in the output to explain it. An
@infoon fallback (suggestion 2) is the minimum; recording the solver in the Galerkin HDF5 group would be better. - Accuracy at this head. Could you post the Cholesky-vs-LU difference at
c827d4d98, plus a comparison that says which is closer to converged (the Riccati Δ′ on the same equilibrium, or agal_nxrefinement)? A 5.6e-5 relative change in the diagonal moved Δ by 0.4 %, so the extraction is sensitive. - Split out
extended_precision_bvp. Double64 → Float64 is a precision change unrelated to the Galerkin solver, and it also changes whatctrl === nothingcallers get. I'd move it to its own PR with isolated harness numbers (diiid_n1_riccati,solovev_*, a small-singfac_mincase), the timing behind "1.5–5.4×", and a test of the imaginary and off-diagonal entries (suggestions 5–7). - Solve-stage timing for Cholesky, LU, and a failed Cholesky followed by LU. On a free-boundary case that is indefinite only in the edge vacuum block, the failed attempt costs close to a full factorization.
- Merge develop and re-run the harness. Develop has since changed the Galerkin element sizing and the GGJ ray path in the same files, and still carries
alpha_override, so the current report is against a stale base.
Smaller: the release note's item 1 quotes numbers from before the M₀-block change; the kl != ku case reports "not positive definite" for what is a structural mismatch; and the new resonant-trace warning fires per side on every crossing and triggers on any small trace near Mercier-marginal surfaces, with an undocumented 1e-3 threshold.
| fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n] | ||
| ab = zeros(ComplexF64, kl + 1, n) | ||
| for j in 1:n, i in j:min(n, j + kl) | ||
| ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j) | ||
| end | ||
| gal_zpbtrf!(ab, kl) || return false |
There was a problem hiding this comment.
Cholesky reads only the lower band and Re of the diagonal, so a non-Hermitian matrix (or a Dirichlet row the pattern match misses) silently solves a different system. Fall back to LU unless the free block is Hermitian to round-off; the negated <= also sends NaN to LU.
| fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n] | |
| ab = zeros(ComplexF64, kl + 1, n) | |
| for j in 1:n, i in j:min(n, j + kl) | |
| ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j) | |
| end | |
| gal_zpbtrf!(ab, kl) || return false | |
| fixed = [all(j == i ? a(i, j) == 1 : iszero(a(i, j)) for j in band(i)) for i in 1:n] | |
| # Cholesky reads only the lower band: fall back unless the free block is Hermitian to round-off. | |
| for j in 1:n, i in (j+1):min(n, j + kl) | |
| fixed[i] || fixed[j] || abs(a(i, j) - conj(a(j, i))) <= 1e-10 * sqrt(abs(a(i, i)) * abs(a(j, j))) || return false | |
| end | |
| ab = zeros(ComplexF64, kl + 1, n) | |
| for j in 1:n, i in j:min(n, j + kl) | |
| ab[1+i-j, j] = fixed[i] || fixed[j] ? (i == j ? one(ComplexF64) : zero(ComplexF64)) : a(i, j) | |
| end | |
| gal_zpbtrf!(ab, kl) || return false |
There was a problem hiding this comment.
Implemented the off-diagonal guard, including finite-entry checks and Dirichlet-row handling, in 633f9bd and retained it in the final head. Entries use the local scale 1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj)); rejected attempts preserve the full matrix/RHS for LU. One deliberate qualification: the diagonal is still projected to its real part, so this is not a full Hermiticity check. That choice is explicit in the docs and a separate imaginary-diagonal test. The coefficient alignment subsequently reduced the observed LU/Cholesky Delta-prime difference to at most 2.53e-9 relative norm across our expanded battery. See implementation and tests.
| function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String) | ||
| if solver != "LU" | ||
| gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky" | ||
| solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not positive definite; use \"auto\" or \"LU\"") | ||
| end | ||
| gal_scaled_lu_solve!(mat, sol, kl, ku) | ||
| return "LU" | ||
| end |
There was a problem hiding this comment.
Report a fallback even without verbose, so a solver switch inside a scan is visible, and make the forced-Cholesky error name both failure modes.
| function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String) | |
| if solver != "LU" | |
| gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky" | |
| solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not positive definite; use \"auto\" or \"LU\"") | |
| end | |
| gal_scaled_lu_solve!(mat, sol, kl, ku) | |
| return "LU" | |
| end | |
| function gal_banded_solve!(mat::Matrix{ComplexF64}, sol::AbstractVecOrMat{ComplexF64}, kl::Int, ku::Int, solver::String) | |
| if solver != "LU" | |
| gal_cholesky_solve!(mat, sol, kl, ku) && return "cholesky" | |
| solver == "cholesky" && error("gal_solver=\"cholesky\": the Galerkin matrix is not Hermitian positive definite; use \"auto\" or \"LU\"") | |
| @info "Galerkin matrix is not Hermitian positive definite; solving by LU" | |
| end | |
| gal_scaled_lu_solve!(mat, sol, kl, ku) | |
| return "LU" | |
| end |
There was a problem hiding this comment.
Implemented: auto fallback emits an info message even with verbose=false, and the forced-Cholesky error names off-diagonal symmetry/real-diagonal positive-definiteness failure. Unequal kl/ku gets its own structural diagnostic. We also went beyond logging: GalerkinResult records the factorization actually used, written to ForceFreeStates/Solutions/GalerkinIntegration/solver. The physical free-boundary timing run verified an actual failed banded Cholesky followed by LU, rather than only rejection by the guard.
| @test_throws ErrorException FFS.gal_banded_solve!(ab, copy(rhs), kl, kl, "cholesky") | ||
| end |
There was a problem hiding this comment.
A matrix whose lower band is positive definite but which is not Hermitian must go to LU. This fails on the current head.
| @test_throws ErrorException FFS.gal_banded_solve!(ab, copy(rhs), kl, kl, "cholesky") | |
| end | |
| @test_throws ErrorException FFS.gal_banded_solve!(ab, copy(rhs), kl, kl, "cholesky") | |
| # Positive definite lower band but non-Hermitian: Cholesky would solve a different system. | |
| H[5, 5] = 4 + sin(5) | |
| H[6, 4] += 0.1 | |
| A, ab = system(H) | |
| sol = copy(rhs) | |
| @test FFS.gal_banded_solve!(ab, sol, kl, kl, "auto") == "LU" | |
| @test sol ≈ A \ rhs rtol = 1e-12 | |
| end |
There was a problem hiding this comment.
Added this regression case: a positive-definite lower band with an asymmetric off-diagonal entry must select LU and agree with the dense solution. We also cover a perturbed Dirichlet row, nonfinite entries, unequal bandwidths, tolerated off-diagonal roundoff, and the explicitly retained imaginary-diagonal projection. See the current solver tests. The focused runs passed 27 solver assertions plus 71 coefficient assertions per example.
| # Lower-band Hermitian Cholesky as rdcon gal.f: zpbtrf('L') then zpbtrs('L'). LinearAlgebra.LAPACK | ||
| # has no pbtrf!/pbtrs! wrapper, so call LAPACK directly. Returns false if not positive definite. |
There was a problem hiding this comment.
New backend comment cites a Fortran file, the repo cites papers only.
| # Lower-band Hermitian Cholesky as rdcon gal.f: zpbtrf('L') then zpbtrs('L'). LinearAlgebra.LAPACK | |
| # has no pbtrf!/pbtrs! wrapper, so call LAPACK directly. Returns false if not positive definite. | |
| # Lower-band Hermitian Cholesky (LAPACK zpbtrf/zpbtrs, 'L'). LinearAlgebra.LAPACK has no | |
| # pbtrf!/pbtrs! wrapper, so call LAPACK directly. Returns false if not positive definite. |
There was a problem hiding this comment.
Applied the suggested backend comment; it now describes the lower-band LAPACK zpbtrf/zpbtrs wrapper and failure behavior without the Fortran-file reference. See the comment.
| intr, fm_propagators, fm_chunks; | ||
| wv=vac.wv, psio=equil.psio, | ||
| S_at_surface_left=fm_S_left, ctrl=ctrl_ext, equil=equil, mats=mats) | ||
| @test isapprox(real.(diag(intr.delta_prime_matrix)), real.(diag(dpm)); rtol=1e-8) |
There was a problem hiding this comment.
The new precision test checks only the real diagonal, the imaginary and off-diagonal Δ′ were the reason for Double64. Compare the whole complex matrix too.
| @test isapprox(real.(diag(intr.delta_prime_matrix)), real.(diag(dpm)); rtol=1e-8) | |
| @test isapprox(real.(diag(intr.delta_prime_matrix)), real.(diag(dpm)); rtol=1e-8) | |
| @test isapprox(intr.delta_prime_matrix, dpm; rtol=1e-8) |
There was a problem hiding this comment.
Implemented the whole-complex-matrix comparison in addition to the real diagonal, and added per-entry real/imaginary checks so a dominant large entry cannot mask every small component. The test explicitly compares the retained default Double64 path with extended_precision_bvp=false. A separate helper test uses an independent 256-bit reference. See the integration comparison and the reference test.
| integrator::String = "riccati" # Default: unlocks SingularSurfaces/Delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. Use "forward" for dense ξ (PerturbedEquilibrium) or kinetic runs. | ||
| nchunks::Int = 0 # Riccati chunk-count target; 0 = auto (derived from msing alone, never from Threads.nthreads()). | ||
| extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). | ||
| extended_precision_bvp::Bool = false # Δ' BVP in Complex{Double64}; opt-in check of Float64 at small singfac_min. |
There was a problem hiding this comment.
Keep the Double64 default here and move the flip to its own PR with isolated numbers (restore the line-156 docstring with it).
| extended_precision_bvp::Bool = false # Δ' BVP in Complex{Double64}; opt-in check of Float64 at small singfac_min. | |
| extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). |
There was a problem hiding this comment.
Done in 633f9bd: extended_precision_bvp=true is restored, together with the default/precision contract in its docstring. This PR does not flip that default; any proposed precision reduction will need separate evidence and review. See the control definition.
| end | ||
| # Fortran sing_solve: vmat(n1,isol,:,k) *= sig/(power(isol)+k/2) | ||
| vmat[n1, isol, :, k+1] .*= sig / (power[isol] + k / 2.0) | ||
| end |
There was a problem hiding this comment.
Whatever the default, a ctrl === nothing caller should keep getting what it got before.
| Tc = (ctrl === nothing || ctrl.extended_precision_bvp) ? Complex{Double64} : ComplexF64 |
There was a problem hiding this comment.
Done: the helper chooses Complex{Double64} when ctrl === nothing or the option is true, preserving the old no-control behavior. We also added a no-control test against a 256-bit reference, with an explicit Float64 comparison. The helper now lives in DeltaPrimeBVP.jl.
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
One structural point on top of today's review, since it is not in the threads: the trace-aware Frobenius exponents (4b44dcd, the former #517) are what moves the numbers here (Im Δ′ at q = 5 from −15182 to 24.9, the 82 % matrix-norm change), they are agreed correct above, and the harness table shows the later review commits change nothing on top of them. The Cholesky items still open (Hermiticity before auto, the timing and accuracy data, the develop re-run) do not touch that commit. Splitting the exponent fix into its own small ForceFreeStates - BUGFIX! - PR with its own harness stamp would land the physics correction now and leave this PR to the solver question; it also unblocks the parity plan in #591 (item A7). At 633f9bd the Hermiticity guard and the visible fallback are in and extended_precision_bvp is back to true, so of the six items only the three data requests and the develop merge remain on the author.
|
@logan-nc fair points but I think we have discussed and agreed on the fortran side of things that the fallback solver model is best, and we are going to go ahead and implement that now. I will merge this in the next couple of hours and both will go in today. |
|
Summary of how the final branch differs from Stuart's original contribution, and what we added during review/integration: Stuart's contribution remains the foundation. The original submission ( On top of Stuart's latest revision (
What the investigation establishes: matching the coefficients reduces the DIII-D full complex PEST-3 Δ′ discrepancy against Riccati from about 10.98% to 0.152% at nx=256. Across the broader resolution/width/order/quadrature/reference-control battery, LU and Cholesky agree within 2.53e-9 relative Δ′ matrix norm. This supports the coefficient correction and solver consistency; it does not establish Riccati as an exact answer. The residual gap at nx=1024 is about 0.1451% for DIII-D and 0.00207% for LAR. The unresolved accuracy investigation is deferred to #603 by maintainer agreement. Fallback timing is now measured: on a physical DIII-D free-boundary assembly, Cholesky actually fails and auto falls back to LU. Seven warmed, alternating samples give 0.898 s for auto versus 0.540 s for direct LU (+0.358 s, 66%); resulting Δ′ matrices are identical. LAR's free-boundary matrix accepts Cholesky and auto takes 0.299 s versus 0.322 s for LU. These exclude assembly and input copies. Final validation: the post-develop-merge fanout ran all 21 cases at develop, the pre-merge PR head and The PR body contains the final harness stamp and retains the historical investigation record. Third-party human approval is required before merging into develop. |
Release note
6d655c74c, all 63 post-merge runs completed successfully (21 cases × develop6d655c74c, pre-merge PRc00ede00e, merged heada1efcd6f7). Compared withc00ede00e, 20 cases are unchanged within thresholds and the driven-PE case retains eight N/A quantities; no quantities are flagged changed by the merge. In the preceding coefficient-change comparison (b75eeb3a8→c00ede00e), 18 cases are unchanged within harness thresholds, two Galerkin cases each have five changed quantities, and the driven-PE case retains eight N/A quantities on all refs. In both Galerkin cases, the coefficient correction changes the PEST-3 Δ′ diagonal by up to 12.16%, the raw Δ′ matrix norm by +5.38%, the coil-block norm by −93.51%, matchedcoutnorm by −0.43%, and matchedbpennorm by −67.06%. These are material driven-response changes requiring reviewer attention, not a claim of independently validated absolute accuracy.gal_solverdefaults to"auto", which attempts Cholesky and falls back to scaled LU. Explicit"cholesky"requires Cholesky to succeed and errors otherwise;"LU"skips the attempt. Selecting LU does not undo the Frobenius or coefficient corrections.extended_precision_bvpretains defaulttrue;falseremains an explicit Float64 option. Re-baseline affected Galerkin Δ′ and driven-response results.The combined change implements banded Cholesky with checked failure and LU fallback, preserves nonzero Dirichlet boundary values, and retains the trace-aware Frobenius exponent correction. Galerkin now reconstructs F̄=LL† from the same interpolated factor L used by the Euler–Lagrange kernel and singular expansions, rather than independently interpolating F̄. Existing buffers are reused for the reconstruction.
Implementation and focused validation
Published head:
a1efcd6f70b0d888b90415de6da191e7af704146. Develop6d655c74cwas merged without conflicts; all 29 merged Julia files parse and whitespace checks pass. The worktree is clean. The full harness was rerun at this merged head.1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj)). It deliberately projects the diagonal onto its real part; this is not a full-matrix Hermiticity test.verbose=false; unequal bandwidths receive a structural diagnostic. The actual factorization is recorded underForceFreeStates/Solutions/GalerkinIntegration/solver.Δ′ comparison and deferred accuracy investigation
The comparison is the full complex PEST-3 matrix norm
norm(Δ′_Gal - Δ′_Ric)/norm(Δ′_Ric), including off-diagonal coupling. Both methods use the same equilibrium, Fourier band, surface ordering and fixed displacement boundaries; Riccati uses the extended-precision BVP. These comparisons do not establish Riccati as an exact reference.Aligning the coefficient representation reduces the DIII-D discrepancy at nx=256 from about 10.98% to 0.152%. Across the expanded battery, LU and Cholesky agree within 2.53e-9 relative Δ′ matrix norm.
An 11-Galerkin/6-Riccati-variant battery per equilibrium varied radial resolution, matching widths, expansion order, quadrature, Riccati tolerance and crossing distance. A separate eight-variant study isolated the resonant cutoff, larger quadrature budgets and spline-knot splitting. Expansion order barely affects the residual, and Riccati-control effects are much smaller than it. Cutoff/budget changes do not explain the gap. All 96 instrumented resonant integrals exhausted their quadrature budgets without meeting estimated error targets; increasing the budget 40× barely changes Δ′. DIII-D's resonant intervals contain no interior spline knots, so splitting there was a no-op.
The remaining discrepancy, quadrature stopping behavior, and ordinary/extension-cell knot-crossing investigation are explicitly deferred to #603 by maintainer agreement. The outermost DIII-D surface dominates the matrix norm. Full numerical agreement has not been demonstrated.
Factorization timing
Physical free-boundary assemblies at nx=256, without RPEC coil columns. Medians of seven warmed samples, alternating solver order; assembly and input copies excluded. Auto includes the guard, Cholesky attempt and any LU fallback.
DIII-D's resulting Δ′ matrices are identical between direct LU and fallback. LAR's relative difference is 7.74e-11. Auto is not universally faster: known indefinite systems can select
"LU"directly.Final post-merge regression harness: a1efcd6
21 cases × three refs: develop
6d655c74c522313084beb3807a8f71a78f28785a, pre-merge PRc00ede00e538596f10ef0fb4e791de4c11d3e6c3, merged heada1efcd6f70b0d888b90415de6da191e7af704146. All 63 case runs, preparation and report completed successfully; report exit code 0. Julia 1.12.6, CI-pinned manifest4c9d54fa, Slurm fanout, matching per-case resources across refs.Compared with the pre-merge PR, 20 cases are unchanged within harness thresholds, including both Galerkin cases.
gal_resistive_peretains eight N/A quantities on all refs. No tracked quantities are flagged changed by the develop merge. The material coefficient/driven-response changes established below remain part of the combined PR; the N/A case remains a coverage limitation.Coefficient-change harness before the develop merge: c00ede0
21 cases × three refs, Julia 1.12.6, CI-pinned manifest
4c9d54fa, executed through Slurm fanout. All tasks and the final report completed successfully; report exit code 0. Per-case resources were identical across refs, including the separate precompiled Riccati case.6d655c74cb75eeb3a8c00ede00ecoutnormbpennormSurface locations, D_I, α and inner-layer Δ remain unchanged by the final coefficient commit. Matched residuals remain near 1.7e-16. The harness flags the five changed quantities in each Galerkin case; an exit code of zero does not mean all numerical results are unchanged. The driven-PE case's eight missing quantities remain a validation limitation.
Review status
The implementation, convergence battery, fallback timing and final harness are complete. Residual accuracy investigation is tracked in #603. The large coil/driven-response changes above remain part of the human review scope.
Third-party human reviewer approval is required before merging into develop.
Historical PR body and investigation reports (superseded by the final results above)
Release note
4c4126738, the previous PR head74843b814, and review-fix head633f9bd33. The combined changes reduce the Galerkin raw Δ′ matrix norm from about 12438 to 2184 and change matchedbpennorms by 3.11%. The review fixes leave tracked Galerkin quantities unchanged against74843b814; Riccati diagonal differences are about 1.6e-10, and SLAYER growth rates change by 0.02%. This stamp predates the latest develop merges. The final harness will run after the remaining accuracy and timing work is finished.gal_solvernow defaults to"auto"(Cholesky first, then LU); set"LU"to select the previous factorization. That selection does not undo the Frobenius exponent correction.extended_precision_bvpretains its previous defaulttrue, including for callers without a control object;falseremains an explicit Float64 option. Results pinned to the old imaginary Δ′ or Galerkinalphaneed re-baselining.The Galerkin solver now uses working banded Cholesky wrappers, checks factorization failure, and falls back to scaled LU. Dirichlet columns are transferred to the RHS, preserving nonzero coil boundary values. The trace-aware Frobenius exponents from #517 remain in this PR: each side uses the eigenvalues of its own resonant M₀ block, with consistent root ordering on Mercier-unstable surfaces.
Current implementation and integration status
b75eeb3a8, incorporating develope0324a91d(including the SLAYER transport-model updates). This merge was clean; all 22 merged Julia files parse and whitespace checks pass. The convergence measurements below remain measurements at3d3e723ef.1e-10 * sqrt(abs(A_ii)) * sqrt(abs(A_jj)). The diagonal is deliberately projected onto its real part; this is not a full-matrix Hermiticity test.verbose=false. Unequal bandwidths get a structural diagnostic. The factorization actually used is written toForceFreeStates/Solutions/GalerkinIntegration/solver.abs(tr(M₀))/2 > max(1e-3*abs(α), 64*eps(Float64)*opnorm(M₀, Inf)). It diagnoses trace consistency, not radial/Fourier convergence.3d3e723ef, GitHub's Julia 1.11 and latest-Julia test suites, aggregate Tests check, and documentation pass. The release-note check remains stale until the final harness stamp is available.Convergence and successful-solve timings at 3d3e723
Completed on ITER SDCC (Slurm array 1296767, both tasks exited 0), Julia 1.12.6, pinned manifest
4c9d54fa, four Julia threads and one BLAS thread per task.Both shipped examples were evaluated at
gal_nx = 128, 256, 512. LU and Cholesky-first used identical assembled matrices and RHSs. The Riccati references used the same equilibrium, Fourier band, surface ordering, and fixed outer displacement boundary condition.At
gal_nx = 256:Solver differences are
norm(Δ′_auto - Δ′_LU)/norm(Δ′_LU); refinement differences use each solver's nx=512 result as the denominator. Timings are medians of three warmed solves, excluding input copies and assembly, and including the guard, factorization, and multi-RHS solve. Cholesky was actually used in all six Galerkin runs. The largest DIII-D-like relative change in a real diagonal entry between solvers at nx=256 is 2.65e-7.Accuracy limitation: the Galerkin–Riccati discrepancy remains about 11% on DIII-D-like and 0.0558% on LAR at nx=512. At DIII-D-like q=5 and nx=256, LU gives approximately −2703.9317 + 25.1123i, Cholesky gives −2703.9324 + 25.0880i, and Riccati gives −3036.3550 + 2.1e-5i. Refining only the Galerkin element count does not remove this difference. These results establish close LU–Cholesky agreement on the tested systems; they do not establish which formalism is closer to the continuum result.
Remaining validation before the final harness
Last complete regression harness: 633f9bd
regress --refs 4c4126738f4a5d1a473b9814d4e661113542698f,74843b814f4677536f09b003b250f00248e5212a,633f9bd33568e23d22a7928f042d71dea5f88287All 57 runs (19 cases at each of three refs) completed successfully on ITER SDCC with Julia 1.12.6 and the pinned CI manifest. The comparison of
633f9bd33with74843b814reports 17 cases unchanged within their harness thresholds, one case with a 0.02% SLAYER growth-rate change, andgal_resistive_pewith all eight quantities N/A on both refs. Tracked Galerkin quantities are unchanged; the Riccati BVP diagonal's maximum absolute difference is about 1.6e-10. The 3.11% matched-bpen shift against develop remains part of the combined PR, not a new shift introduced by the review fixes.Historical reports
The earlier reports below are retained as the investigation record. Their heads, defaults, and validation statuses describe those earlier revisions.
Earlier harness reports and measurements
Regression report
regress --refs 23ff553a5140d813f5d76549b70b3fcaa74b1769,74843b814completed on ITER SDCC with Julia 1.12.6 and the pinned CI manifest4c9d54fa. All 38 runs (19 cases at each ref) completed successfully. Most cases used 2 Julia threads / 1 BLAS thread;diiid_error_fieldused 16 Julia threads / 8 BLAS threads on both refs.gal_resistive_diiid,gal_resistive_diiid_galinnerdiiid_n1_riccati,diiid_n1_riccati_precompileddiiid_n1diiid_multi_ndiiid_slayer_n1,efit_fixedbdy_separatrix,ggj_ray_q500i,ggj_reference,solovev_kinetic_calculated,solovev_kinetic_nuzerogal_resistive_peThe full Julia test suite at
74843b814passed on ITER (Julia 1.12.6, 4 CPUs, Slurm job 1296033), with one test marked broken. GitHub's latest-Julia tests and documentation also passed. The cancelled Julia 1.11 job has been resubmitted.Historical regression report
regress --refs 278057d96,c827d4d98(develop vs this head, all 18 cases; julia 1.11.7, pinned manifest, 2 threads):gal_resistive_diiiddiiid_n1_riccati(and_precompiled)diiid_n1,diiid_multi_n,diiid_truncation_n1,diiid_error_fieldsolovev_n1,solovev_multi_n,solovev_kinetic_*diiid_slayer_n1gal_resistive_peis N/A on both, as on developWith
regress --refs 278057d96,4b44dcd02,c827d4d98, the three review commits (a968945, 2c02244, c827d4d) change none of these numbers: every case reports c827d4d unchanged against 4b44dcd.Galerkin Cholesky alone, before #517 was merged (6fc89d0 vs cd290c6)
regress --cases gal_resistive_diiid,gal_resistive_pe --refs 6fc89d021,cd290c659(develop merge-base vs the Cholesky commit; julia 1.11.7, pinned manifest, 2 threads).gal_resistive_peis N/A on both refs, as on develop.